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ABSTRACT 

We present a set of global, self-consistent iV-body/SPH simulations of the dynamic evo- 
lution of galactic discs with gas and including magnetic fields. We have implemented a 
description to follow the evolution of magnetic fields with the ideal induction equation 
in the SPH part of the Vine code. Results from a direct implementation of the field 
equations are compared to a representation by Euler potentials, which pose a V-B-free 
description, an constraint not fulfilled for the direct implementation. All simulations 
are compared to an implementation of magnetic fields in the Gadget code which 
includes also cleaning methods for V • B. 

Starting with a homogeneous seed field we find that by differential rotation and 
spiral structure formation of the disc the field is amplified by one order of magnitude 
within five rotation periods of the disc. The amplification is stronger for higher numer- 
ical resolution. Moreover, we find a tight connection of the magnetic field structure 
to the density pattern of the galaxy in our simulations, with the magnetic field lines 
being aligned with the developing spiral pattern of the gas. Our simulations clearly 
show the importance of non-axisymmetry for the evolution of the magnetic field. 

Key words: methods: iV-body simulations - galaxies: spiral - galaxies: evolution - 
galaxies: magnetic fields - galaxies: kinematics and dynamics 



1 INTRODUCTION 

Radio observations have revealed that disc galaxies are per- 
meated by large scale mag netic field s ordered on kpc scales 
and beyo nd (Be ck &; Hoe rnes 19961 iHummel Beck 1995, 
iBeck et al.l [1985). The typical field strength, determined 
from polarization, Faraday rotation and energy equiparti- 
tion is of the order of 10 /xG (e.g. iBeckl 120041 ). The spa- 
tial structure of the B-field reflects the spiral and/or barred 
structure o f the gas distribution within the galactic discs 
(|Beckll2008l V For example, Fig. [T] shows optical observations 
of the spiral galaxy M51 overlayed with contours of total 
synchrotron intensity (tracing the total magnetic fleld) and 
magnetic field vectors. It reveals the tight connection of mag- 
netic fleld with the gas distribution in the galactic disc. 

The motion of the gas within the gravitational potential 
of a galaxy strongly influences the strength and direction 
of the magnetic fleld in the interstellar medium. This can 
be seen by inspecting the well known induction equation of 
magnetohydrodynamics (MHD), i.e. the temporal evolution 
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equation for the magnetic fleld, 

^ = V X (v X B) - (V X ry(V X B)), (1) 

where v denotes the gas velocity and r] represents the mag- 
netic diffusivity which is inversely proportional to the elec- 
trical conductivity. 

Apparently, within the frame of MHD, the role of the 
galaxy as a whole is simply to provide for the gas velocity 
fleld. Since the conductivity of the interstellar medium is 
very high, the magnetic fleld is closely coupled to the gas 
motion. It is this 'frozen-in'-property of both, the magnetic 
fleld and the gas, which determines the spatial structure of 
the magnetic fleld. In other words, a detailed investigation 
of the velocity fleld of the interstellar gas in disc galaxies is 
necessary for a deeper physical understanding of the evolu- 
tion of galactic magnetic flelds. 

The gas in the disc rotates differentially within the 
global gravitational potential. Angular momentum trans- 
port via spiral arms, bars and gravitational interaction forces 
the gas to move towards the central regions, and eventually, 
star formation activity in the disc (superbubbles, winds etc.) 
drives gas perpendicular to the plane of the disc towards the 
galactic halo. In general, the axisymmetric rotation velocity 
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Figure 1. Optical image of M51 (Hubble) overlayed with con- 
tours of total synchrotron intensity as measure for the total 
magnetic field (combined observations at Effelsberg and VLA 
at 6 cm) and vectors of magnetic field. From A. Fletcher &; R. 
Beck (MPIfR) and Hubble Heritage Team (STScI), published by 
'Sterne und Weltraum', September 2006. 



is the dominant component, followed by non-axisymmetric 
and radial components. The velocity components perpen- 
dicular to the disc are typically the smallest. Altogether, 
V in Eq. [T] represents a complex three-dimensional non- 
axisymmetric velocity field strongly coupled to the global 
properties of the galaxy, including the dark matter halo, 
stellar component and internal disc activity. 

Beside the large scale components of the gas veloc- 
ity field there are also small scale velocity fluctuations of 
interstellar gas driven by all kinds of local disc activity, 
i.e. stellar winds, supernova expl osions, cloud-cloud colli- 
sions , galactic winds, etc (see e . g. lFerrier3 19921, Efstathiou 
200C, Johansson Efstathioul l2006l . I Kulsrud Zweibell 
2008 , iGressel et al.l l2008l ) . These unordered velocity com- 
ponents generate two effects which are known as helicity 
(in terms of a convective turbulent motion perpendicular 
to the disc) and turbulent diffusion (magnetic field lines 
with antiparallel direction reconnect and annihilate par- 
tially). Helicity supports the amplification of the magnetic 
field, whereas turbulent diffusion reduces the field strength 
(see, e.g. Brandenburg & Subramanian 2005 for a review 
of nonlinear dynamo theory). Therefore, an incorporation 
of these small scale velocity components into the analy- 
sis requires some manipulation of th e induction equation 
(Eq. [1]) in terms of a mean- field theory feteenbeck &: Krausd 



Il969l . IWielebinski k KraT[s3ll993l . [s"ur et al.|[20Q7h . Within 
the frame of the mean-field description the velocity and mag- 
netic fields are considered as superpositions of the mean 
and fluctuating parts (v = (v) + and B = (B) + B^). 
The fluctuating velocity components are coupled to small- 
scale fluctuations of the magnetic field. The coupling terms 
are then given by V x a(B), where a = |^(v^ • (V x v^)) 
f Zeldovich et al.l ll983\ and by ryTA(B), where r/T now de- 
scribes the turbulent diffusion coefficient r/r oc t'turb • ^turb, 
where Vtuvh and /turb are the typical velocity and length scale 
of the turbulent motion, respectively. 
This leads to the dynamo equation 



^B 

— = Vx(vxB)+Vx aB, 



(2) 



where we have neglected the diffusivity and dropped the 
mean-brackets for convenience (here, and in the following, 
B and v refer to their mean values). 

Eq. [2] is the central equation of cosmic mean field dy- 
namos. It describes the circle of amplification of the dif- 
ferent components. The classical dynamo model describes 
the amplification of the magnetic field through the follow- 
ing chain of a (convective turbulence) and Q (differential 
rotation) actions: The radial component Br is amplified 
through a-action from turbulence; then is generated 
from Br through Q-action from the shear of the galactic 
differential rotation. Such an aQ mean field dynamo am- 
plifies the mag netic field by repe ating the chain of a and 
Q actions (see IWidrowl I2OO2I and IStefani et al. I I2OO8I for a 
review of dynamo theory). However, the origin of the a- 
effect is still under discussion dCattaneo Vai nshtein 199lL 
IVainshtein Cattaneolll99llKulsrud Ande rson 1992^). 

We emphasize that the described classical dynamo mod- 
els use only one velocity component, the differential rotation. 
To be more precise the role of any deviation from axisym- 
metry is considered to be unimportant for the evolution of 
the large-scale magnetic field, which is not necessarily true 
in real galaxies. 

On this account, there have been three-dimensional 
numerical simulations using an analytical turbulent veloc- 
ity field, where deviations from axisym metry were incorpo- 
rated in the gas- and t urbulence-profiles (|Rohde et al.ll 19971 , 
iRohde Elstnerlll998l ). These studies showed, that even ac- 
counting for the a-effect calculated out of the analytical ve- 
locity field an initial magnetic field cannot survive for more 
than 500 Myr. 

Moreover, lElstner et al.l I2OOOI performed iV-body sim- 
ulations of two component (collisionless stars and gaseous 
clouds moving in the gravitational potential of the stellar 
population) , self- gravitating discs embedded in an analytical 
bulge- and halo-potential. These simulated clouds provided 
an already very good approximation of the gas velocity field. 
However, full hydrodynamics was not incorporated. The ob- 
tained velocity field was used in an aQ-dynamo description. 
Without including the a-effect, the non-azimuthal 3D gas 
flow alone did not provide an ampliflcation of the magnetic 
fleld. The fleld got amplifled by several orders of magnitude 
within 0.7 Myr only when the a-effect was included. In ad- 
dition, they found an alignment of the magnetic fleld with 
the developed spiral pattern o f the d isc. 

Recently, iDobbs Pried (|2008l l performed three di- 
mensional, full MHD, single and two-component (cold and 
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hot gas) simulations using smooth particle hydrodynamic 
(SPH) methods to treat MHD. They applied a spiral po- 
tential to the gas, thus, the self-induced formation of spi- 
ral structure was not included. Their work concentrated on 
structure formation in the disc, like molecular clouds and 
inter- arm spurs. They found that the main effect of adding 
a magnetic field to these calculations was to inhibit the for- 
mation of structure in the disc. They did not consider global 
enhancement and structure formation of the magnetic field, 
but nevertheless, they found that the global magnetic field 
was following the large scale velocity field. 

It is the aim of this paper to present further steps to- 
wards a more complete dynamo model. We perform for the 
first time a set of self-consistent iV-body calculations of a 
spiral galaxy including hydrodynamics as well as the induc- 
tion equation via the SPH method to obtain the complex 
three dimensional velocity field. Compared to all previous 
work, we use no analytical potential for any component of 
the galaxy. All components (disc, gas, bulge and halo) are 
represented by particles which are treated as self- gravitating 
AT-body-particles, while hydrodynamics is applied to the gas 
component only. We use more than o ne order of magni- 
tude more particles than lElstner et al] (|2000h . We follow 
the evolution of the magnetic field according to the induc- 
tion equation (eq. [T]). Thus, we have implemented the SPH 
variant of the induction equation as well as the representa- 
tion of the magnetic fields by Euler potentials in the SPH 
code Vine and compare the results with simulations per- 
formed using the SPH code Gadget. AT-body/SPH methods 
are well adapted for simulating whole galactic discs as the 
simulated discs stay stable for at least 15 dynamical times 
(where we define the dynamical time for a disc galaxy as its 
half mass rotation period). As we will show in section [4] our 
discs are forming spiral structure without applying a spiral 
potential or any other mechanism to provide extraordinary 
flows. 

In summary, we investigate here the kinematic reac- 
tion of a large-scale magnetic field on the complete three- 
dimensional, large-scale velocity field of a disc galaxy ob- 
tained from the AT-body SPH simulations, using two differ- 
ent numerical codes. 

The paper is organized as follows: Sect ion [2] gives shortly 
the theory of magnetic field evolution in differentially ro- 
tating systems. A summary of the SPH method and the 
treatment of magnetic fields including the method based on 
Euler potentials is given in section [S] The simulations to- 
gether with a comparison of the performance of the Vine 
and Gadget codes are presented in section |4l The results 
are discussed in section \5\ where we also analyse the terms 
of the induction equation in detail. Finally, we summarise 
and conclude in section (6] 



2 THEORETICAL EXPECTATIONS 

When only studying the effect of the gas velocity on the evo- 
lution of the magnetic field, we can neglect the diffusive term 
in eq. [T] Keeping this term, one would physically except the 
magnetic field to dissipate depending on the value of 77 and 
reconnect when converse magnetic field lines come together. 
Technically, 77 is not always assumed to be spatially depen- 
dent, so that the diffusive term reads —r]V^'B. However, this 



formulation leads only to an effective smoothing, and not a 
real diffusion of the magnetic field. Neglecting the diffusive 
term thus corresponds to considering an upper limit of field 
amplification. Additionally, 77 is assumed to be small except 
within strong shocks. 

The induction equation [1] then yields 



dB 

dt 



(B • V)v + V (V • B) -B(V • v) - (v V)B, (3) 



and applying cylindrical coordinates eq. [3] reads: 
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These equations can be simplified to get a first idea of how 
magnetic fields will evolve in a galactic disc. For a differen- 
tially rotating disc {dv^/dr ~ 0) with a perfectly axisym- 
metric velocity field, v does not depend on ip. The same 
holds for an axisymmetric magnetic field. If we also assume 
that changes of all quantities in the z direction are small 
compared with those in the radial direction and Bz ^ 0, the 
equations for the regular field, i.e. the field in the plane of 
the disc, read: 



dBr 
dt 

mp 

dt 



-Br - 



-B, 



r 

dvr 
dr 



V^Br o dVr , „ dQ 

— — = -B^^ + Brr—-, 
r or or 



(7) 

(8) 



where Q is the angular velocity. 

In the absence of radial flows, the last term of eq. [8] 
describes the generation of a toroidal magnetic field from 
the radial component of the already present magnetic field 
by differential rotation. This effect is the so called Q-effect 
already mentioned above. Since ^ {vr,Vz) this term is 
dominant and one would expect any initial magnetic field 
to be first wound up by differential rotation. However, this 
effect alone cannot be responsible for a significant amplifi- 
cation of the magnetic field, as the amplification stops when 
all of the radial field is wound up. However, if a gas flow in 
the negative radial direction (i.e. towards the centre of the 
disc) is present the radial field can be amplified and then be 
converted into a toroidal field. These radial gas flows occur 
when angular momentum is transported in the gas out of 
the disc, e.g. by spiral arms or bars. The toroidal magnetic 
field can be amplified further if the radial gas flow velocity 
decreases with increasing galactocentric radius. 

Therefore, a good understanding of the evolution of 
magnetic fields in galactic discs requires full information of 
the three-dimensional velocity field of the gas which is natu- 
rally provided by self-consistent numerical simulations. We 
will discuss the velocity field in our simulations and the re- 
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suiting values of the different terms of the induction equation 
in section [5] 



first law of thermodynamics. This leads to the following SPH 



variant of its ideal form (- 



VP 

p 



3 NUMERICAL METHODS 
3.1 Vine 



The equations presented below are implemented within the 
OpenMP parallel iV-body/SPH evolut ion code Vine . For 
all details we refer t he reader to .Wetzstein et al.l (|2008l ) and 
iNelson et al] (|2008l ). 



3.1.1 SPH basics 

Within the SPH formulation a hydrodynamic quantity A is 
interpolated by a kernel function W{r — r\h) with f Wdr = 
1 and lim/t^o W = S(r — r), where the so called smoothing 
length h defines the spatial extent of the function W. This 
interpolation is then discretised, so that 



Pi 



(9) 



where i (j) is the index of the particle at position (r^) 
and Ai (Aj) the value of the quantity A at the position of 
particle i (j). pj and rrij denote the density and mass at 
position of particle j, respectively. 

The Vine code use s the c ommon W4 kernel defined by 
iMonaghan &; Lattanzid (Il985l) as 



Wij = W{r^j,h) 



else 



(10) 



where values with index ij denote differences (e.g. Yij = 
— Yj) and arithmetic means (e.g. hij — 0.5 • {hi + hj))^ 
respectively, q — |r — r^|//i, v is the number of spatial di- 
mensions of the system and cr is a constan t of order un ity. 
See Monaghan ( 1992 ). .Monaghan (,2QQli ) or lPricd (|2005l ) for 
more details. 



3.1.2 Continuity equation 

As long as the kernel itself is differentiable, every function 
A can be interpolated to a differentiable function by the 
procedure described above. The most common formulation 
of derivatives in SPH is (see e.g. Monaghan 1 992„ . .Price 2005 ) 



{VA)i = — ^mj{Aj - Ai)\/iWi:^ 



(11) 



Using the continuity equation, the total time derivative of 
the density thus reads 



dt 



(12) 



3.1.3 Momentum and Energy equation 

A natural ansatz to derive a conservative form of the mo- 
mentum equation comprising the force due to pressure gra- 
dients (in addition to the force due to the gravitational po- 
tential) is to use the Lagrange formalism together with the 



dt - p. - ^"^'U^ 



(13) 



In this formulation momentum is conserved exactly, since 
the contribution of particle j to the momentum of particle i 
is equal and negative to the contribution of particle i to the 
momentum of particle j. 

The change in the thermodynamical state of the gas 
requires an evolution equation for a state variable corre- 
sponding to the internal energy or entropy of the gas. Vine 
employs an equation for the specific internal energy (u) of 
the gas. Without external heating or cooling terms, only 
compressional heating and cooling are important and the 
SPH variant of the ideal form (^ = — -^V • v) reads: 



duj _ Pi_ ^ 



(14) 



To close the set of equations, an isothermal equation of state 
is used throughout this paper. 



3.1.4 Artificial viscosity 

Artificial viscosity is required to model shocks and angular 
momentum transport properly, where the latter is important 
to be able to simulate spiral arm formation. The Vine code 
uses the most common form of th e artifici al viscosity. It is 
described by the tensor 11^^ as in iMonaghani ((1992 ). 

Since the value of 11^^ depends on the difference in ve- 
locity between the considered particles (i.e. the velocity gra- 
dient) the viscosity increases with increasing velocity gradi- 
ent. Moreover, the viscosity is only applied if particles are 
approaching each o ther. 

iBalsaral (|l995l ) suggested a viscosity limiter to avoid 
spurious angular momentum and vorticity transport in gas 
disks. However, a lower viscosity leads to a higher velocity 
dispersion of the gas and therefore to higher divergence of 
the velocity and magnetic field. As will be shown and dis- 
cussed in section \5\ this higher divergence causes a more 
violent magnetic field amplification. 

The viscous terms within the momentum and energy 
equations read: 



dvi 
~dt 

dui 
~dt 



- - ^ ^ j ^ij V^ Wij 

j 



(15) 
(16) 



This treatment of viscous forces allows for a sensible 
description of the behaviour of gas in a spiral galaxy. How- 
ever, eas. 1141 and 1161 are not applied when using isothermal 
equation of state. 
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3.1.5 Induction equation 



where 



In order to follow the evolution of the magnetic field we have 
additionally implemented equation [3] discretised as 

dt 

^ =-B(V-v) =(B-V)v 



pi . ^ . ^ 



pi 



with /i, V denoting the spatial directions. 



(17) 



3.1.6 Euler potentials 

A well known problem related to magnetic fields within SPH 
is the maintenance of V • B = throughout the simulation. 
Different attempt sjbo s olve this problem have been made 
(see|Price & Monaghan 2005), examples include source term 
approaches ( Powell et al. 1999) and projection methods. 

Theoretically, the problem can be avoided if the mag- 
netic field is represented by Euler potent ials (|Sternl Il97d 
IPrice Bat e 2007, Rossw og Pricell2007l ) . an approach we 
have also implemented into the Vine code. 

The magnetic field is expressed as a function of two 
scalar potentials aE and f3E as: 



B = VaE X VI3e. 

Taking the divergence of B we get: 



(18) 



V-B = V- (Va^; X V/3^;) 
= (V X VaE) - 



-VaE (V X V^e) = 0. (19) 



Thus, the divergence constraint is fulfilled by construction. 

Moreover, for the ideal case (77 = 0) the Euler potentials 
for each particle (i.e. the convective values of the potentials) 
are direct tracers of the magnetic field and stay constant 
with time, 

thus one does not need to perform an additional integration 
when following magnetic fields, leading to a higher accuracy 
of the calculation. The variation of magnetic field is only 
due to the motion of the particles, which corresponds to 
the advection of magnetic field lines by Lagrangian particles 
(frozen flow) (|Sternlll97Ql ). 

Moreover, the formulation by Euler potentials guaran- 
tees conservation of magnetic helicity. 



H 



A Be 



(22) 



which is again reasonable for ideal treatment. Actually, the 
magnetic helicity is zero, since for A = aE^pE and equiv- 
alently A = —[Se^o^e, respectively, B = VaE x Vf3E = 
V X A and therefore A • B = 0. 

The gradients of a^; and Pe can be expressed as 



x7 



'(Vas)r = -;^m,(as,.-a^;,,)(V,m,(/^.))^ (23) 



xr(v/3s)r 



-J2^Ji|3E,^ - |3E,j){V^WiJ{hi)y , (24) 



(25) 



which is exact for linear functions (Price &: Bate|[2007h , i.e. 
for initial conditions with an uniform field. However, the 
difference in using this exact-linear interpolation compared 
to the usual gradient operator is marginal (D. Price, private 
communication) . 

Unfortunately, not all magnetic field configurations can 
be expressed in terms of Euler potentials easily as they enter 
eq.[T8]in a nonlinear way, and they are not unique for certain 
field configurations ("Stern' 'l970l, lYahalom Lvnden-Belll 
I2OO6I ). The former problem restricts only the choice of the 
initial magnetic field, whereas the latter can be crucial. If 
there are field configurations which can be expressed by dif- 
ferent sets of Euler potentials, then this implies, that some 
other field configurations cannot be expressed at all using 
Euler potentials. However, since the fields considered in this 
work are topologically 'simple', we do not expect to en- 
counter these problems. 

Furthermore, Euler potentials do not allow to follow the 
winding of magnetic fields beyond a certain point. This con- 
straint is due to the fact that using the Euler potentials, the 
magnetic field is essentially mapped on the initial particle 
arrangement. If the initial arrangement evolves too much 
during the simulation, particles carrying conflicting values 
of Euler potentials (i.e. values, which do no longer allow 
for a finite and unambiguous calculation of their gradients) 
can come close. Then, the ability of the Euler potentials to 
represent the magnetic field correctly is lost. This conflict 
is expected to occur when the magnetic field is wound up 
more than once, which poses a problem especially towards 
the central region of a simulated galaxy. 



3.1.7 Timestepping 

In Vine, there are basically three different time step criteria, 
based on changes in the acceleration of a particle. 



At: 



its velocity. 



n+l 



At 



or both in combination, 



(26) 



(27) 



(28) 



where e, a and v are the gravitational softening length, the 
acceleration and velocity of a particle in the previous time 
step (n), respectively, and race is an accuracy parameter. 
Two additional time step criteria are applied in SPH simu- 
lations: Fi rst, the Courant- Friedrichs-Lewy criterion as sug- 
gested bv lMonaghanI (|l989l ). 



^^CFL = ^ CFL --TTw 7~75~T V ' (^9) 

Cs -\- l.2{aiCs -\- Pihimaxjfiij) 

where ai and f3i are artificial viscosity parameters, Cs is the 
sound speed, hi the SPH softening length for gas particle i, 
and /Jij corresponds to the velocity divergence between par- 
ticles i and j with the maximum taken over all neighboring 
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particles j of particle i fsee IWetzstein et al.ll2008l for more 
details). Secondly, there is a limit on how much the SPH 
softening lengths are allowed to change during one timestep: 

hi 



hi 



(30) 



where th is again an accuracy parameter. Usually, we apply 
Tacc = 1, TcFL =0.5 and Th = 0.15. The timestep actually 
employed in the simulation is the minimum of the timesteps 
in eqs. [26ll30l 

3.2 Gadget 

A somewhat different treatment of hydrodynamics and 
magnetic fields is realised within the MPI parallel N- 
body/SPH code Gadget (ISpringel et al.1 I2OOI '. 'Springell 
I2OO5I . iDolag k Stasvszvnl I2OO8I 1. There are two significant 
differences in the implementation relevant even for non- 
radiative simulations: 

First, Vine follows a classical implementation which is 
integrating the internal energy, whereas Gadget utilises 
what is generally called the entropy conserving formula- 
tion. The important difference thereby is not the fact that 
Gadget integrates the entropy instead of the internal en- 
ergy. The crucial differences are rather the way in which the 
smoothing length hi is defined (in Gadget, hi is defined 
based on the mass within the kernel instead of the number 
of particles) and the inclusion of correction terms arising 
from the varying smoothing length. Also, the entropy con- 
serving formulation uses a way of symmetrizing the kernel 
given by the derivation of the SPH equations, which in sum 
leads to conservation of ener gy and entropy at the same time 
(|Springel k Hernauistll2002l ). 

The second difference originates in an alternative for- 
mulation of the artificial viscosity. In Gadget, artificial vis- 
cosity is based on the signal velocity instead of sound speed 
(Monaghan 19 97) gind apt to incorporate magnetic waves in 
a natural way ([Price Monaghanll2004 ). 



This different implementation was shown to bring 
measurable improve ments specially for MHD applications 
(jPo lag Stasvszvn 2008), but should not make too much 
of a difference for passive magnetic fields. The implemen- 
tation of the induction equation and the Euler potentials 
formalism is the same in both codes. 

The integration in Gadget is also performed using 
the leapfrog integration scheme, but Gadget utilises a 
kick-drift-kick-scheme whereas Vine uses a drift-kick-drift- 
scheme. 

The timestep is given by 



At" 




(31) 



where 77 translates to the accuracy parameter 

Tacc in eq. 

26l via Tacc = For SPH particles, also a Courant-like 

condition in the form 

Oou^ (32) 



At 



n+1 



is applied, where hi is the SPH softening length for gas par- 
ticle i and v f^ the signal velocity betw een particles i and j 
as defined in lPrice &: MonaghanI (|2004[ ) with the maximum 
taken over all neighboring particles j of particle i. Ccour 



total mass 


Mtot 


= 1.34 • 10 M0 


disc mass 


Mdisc 


= 0.041 • Mtot 


bulge mass 


^bulge 


= 0.01367 • Mtot 


mass of the extended gas disc 


Mgas 


— U.2 • Mdisc 


exponential disc scale length 


Id 


= 3.5 kpc 


scale height of the disc 


h 


= 0.2 -Id 


bulge scale length 


Ib 


= 0.2 -Id 


extent of flat gas disc 


Ig 


6 -Id 


spin parameter 


A 


= 0.033 


virial velocity of the halo 




= 160 km s-i 


half mass circular velocity 


'i^half 


^ 200 km s-i 


half mass rotation period 


^half 


^ 150 Myr 


isothermal sound speed 


Cs 


^ 15 km s"-*^ 


initial magnetic field 


Bo 


10-9 G 



Table 1. Parameters of initial disc setup 



is an accuracy parameter which does not translate one-to- 
one to TcFL in eq. [29] due to the different definition of the 
Courant criterion. We commonly use values of 77 = 0.02 and 
Ccour = 0.15 to ensure that the timestep At in Gadget does 
not get too large compared to Vine. However, changing the 
accuracy parameters by a factor of two does not affect the 
overall evolution and amplification of the magnetic field in 
the simulated systems (not shown). 

Beside that, the codes differ in details of the tree 
construction for calculating gravitational forces. For more 



Wetzstein et a 


[. I2OO8I. iNelson et al. 2008.) 


SDringelll2005l. 


Dolag & Stasvszvn I2OO8I). 



and Gadget 



4 SIMULATIONS 
4.1 Setup 

The initial conditions for our Milky Way like galaxy are re - 
alised using the method described by ISpringel et all (l2005ll 
which is based on lHernquistI (|l993h fsee also I Johansson et al.l 
2009). The galaxy consists of an exponential stellar disc and 
a flat extended gas disc, a stellar bulge and a dark matter 
halo of collisionless particles. The gas is represented by SPH 
particles adopting an isothermal equation of state with a 
fixed sound speed of Cs ~ 15 km s~^, which corresponds 
to a temperature of T 2 • lO'^ K for a molecular weight of 
1.4/1.1 -mproton • We briefly note that by using an isothermal 
equation of state only one component of the ISM is modeled, 
typically this is a reasonably good approximation for the 
warm ga s phase i n disc galaxies (e.g. Barnes 2002, ^ Li et al.l 
I2OO5I . iNaab et al.| [2006 ) . Assuming an isothermal equation of 
state implies that additional heat created in shocks by adi- 
abatic compression and feedback processes (e.g. by SNH) is 
radiated away immediately. In addition, substantial heating 
processes prevent the gas from cooling below its effective 
temperature predefined by its sound speed. 

The parameters describing the initial conditions can be 
found in Table [T] The particle numbers and the gravitational 
and SPH softening lengths used in the different runs can be 
found in Table [2] 

Before we include magnetic fields we allow the galaxy to 
evolve for approximately three half mass rotation periods. 
For simplicity we choose an initial magnetic field in the x 
direction. Its value, Bq = 10 ~^ G, corresponds to the typical 
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low resolution normal resolution high resolution 

number of particles 



Halo 


6 


10^ 


6 • 10^ 


6 • 


10^ 


Disc 


3 


10^ 


3 • 10^ 


3 • 


10^ 


Bulge 


1 


10^ 


1 • 10^ 


1 • 


10^ 


Gas 


3 


10^ 


3 • 10^ 


3 . 


10^ 


Total 


13 


• 10^ 


13 • 10^ 


13 


10^ 






fixed gravitational softening length 


3 e [kpc] 






Vine 


Gadget 


Vine Gadget 


Vine 


Gadget 


Halo 


0.934/2 


0.934 


0.434/2 0.434 




0.199 


Disc 


0.248/2 


0.248 


0.114/2 0.114 




0.052 


Bulge 


0.269/2 


0.269 


0.127/2 0.127 




0.059 


Gas 


0.248/2 


0.248 


0.114/2 0.114 




0.052 






minimum 


SPH softening lengths 


^min 




Gas 


O.Ole 


O.Ole 


O.Ole O.Ole 




O.Ole 



Table 2. Particles numbers and softening lengths. The factor two accounts for the different definition of the Kernel extent in Vine 
{g < 2) and Gadget (g < 1). 
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Figure 2. Surface densities Sgas of the extended gas discs as a 
function of radius before the inclusion of the magnetic fields after 
0.5 Gyr (red line) and after 2 Gyr (black lines) for simulations 
with Gadget (solid line) and Vine (dotted line). The gas discs 
are stable for more than ten half mass rotation periods. 




Figure 3. Same as Fig. [2] but for the stellar surface densities 
Sstars- Both the stellar and the gas discs are stable for more than 
ten half mass rotation periods. 



value of intergalactic magnetic fields feronb erg et al.|[l999[ ) . 
To set up the corresponding Euler potentials, we choose 

ttE = Bo ■ y, (33) 
/3e = y + z. (34) 

We have checked the stability of our discs in indepen- 
dent simulations without magnetic fields. Figs. [2] and [3] show 
the surface densities Sgas of the extended gaseous discs and 
Sstars of the exponential stellar discs, respectively, as a func- 
tion of radius for t = 0.5 Gyr (red), i.e. the time at which the 
magnetic field is switched on, and t = 2.0 Gyr (black). Fig. 
[H shows the circular velocity curves of the simulated galaxies 
at the same times. The discs simulated with Vine (dotted 
Une) and Gadget (soUd line) show similar results and stay 
stable over more than ten half mass rotation periods. 

4.2 Direct magnetic field simulations 

Figs. [5] and [6] show the face on view of the magnetic field 
energy and gas density of the simulated galaxy at different 



output times. The magnetic field was switched on at t = 510 
Myr. The viscosity limiter was not applied. Fig.[5]shows the 
simulation performed with Vine and Fig. [6] the same ini- 
tial conditions simulated with Gadget. The magnetic field 
energy B^/Stt is colour coded and normalised to the initial 
value of ^ • 10~^^ erg cm~^ on a logarithmic scale from 1 
(blue) to 1.5 • 10^ (red). The contours overplotted indicate 
physical densities of 23, 37 and 52 M© pc~^, respectively. 
We use a grid with a cell size of 0.3 kpc for the calculation 
of the mean values of the densities and the magnetic field 
energies, averaging in the vertical direction from -h to /i, 
where h is the local height of the gas disc. 

In both simulations we see that the magnetic field en- 
ergy pattern is tightly connected to the density pattern of 
the gas. Moreover, both simulated galaxies show a very sim- 
ilar morphology in the gas and magnetic field distributions. 
The magnetic field energy in the spiral arms is amplified 
by up to five orders of magnitude in both codes and even 
more in the central region (see also Fig. I13|) . Furthermore, 
the SPH smoothing lengths /igas are similar for both codes 
(Fig. O, indicating that the performance of the hydrody- 
namic calculations is concerted. The smoothing lengths in 
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Figure 5. Face-on magnetic field energy and gas density as a function of time for the simulation performed with Vine using direct 
magnetic field description and without applying the viscosity limiter. The colours correspond to the magnetic field energy B'^ /Stt on a 
logarithmic scale, normalised to the initial value of ^ • lO"-*^^ erg cm~^. The contour lines indicate physical densities of 23, 37 and 52 
Mq pc~^, respectively. 




Figure 6. Same as Fig. [5]for identical initial conditions simulated with Gadget. The morphology is very similar but the magnetic field 
reaches higher values in the spiral arms in the Gadget simulation compared to the simulation with Vine. 



Magnetic fields in spiral galaxies 9 




-ID □ ID -ID ID -1C- a 10 

!t tHpc] K O^SC] M [use) 



Figure 7. Same as Fig. [S] this time the magnetic field is followed using Euler potentials implemented in Vine. In contrast to the direct 
simulation the magnetic field is more strongly amplified in the spiral arms than at the centre. The maximum amplification of the magnetic 
field energy is only three orders of magnitude. Note that the colour scaling is different to Figs. [5] and [6] 
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Figure 8. Same as Fig. [G] this time the magnetic field is followed using Euler potentials implemented in Gadget. The energies and 
morphology of the magnetic field is now similar to the Vine simulation. 
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Figure 4. Circular velocity curves of the simulated galaxies at 
two different times. The colour coding is the same as in Figs. [2] 
and [3] Again, the circular velocity curves are stable over more 
than ten half mass rotation periods. 




Figure 9. SPH softening lengths 
shortly after the inclusion of the magnetic fields at 0.55 Gyr (red 
line), at 1.25 Gyr (green lines) and at 2 Gyr (black lines) for 
simulations with Gadget (solid lines) and Vine (dotted lines). 
The SPH smoothing lengths are very similar for both codes. 



Vine are initially set to a constant value of /igas ~ 0.3 kpc 
at the time of the magnetic field inclusion. 



4.3 SPH with Euler Potentials 

Figs. [7] and [8] show simulations starting from the same initial 
conditions as before. However, this time the evolution of 
the magnetic field was followed using the Euler potentials. 
Again, we show magnetic field energies and gas densities. 
This time the amplification of the magnetic field energy in 
the spiral arms is only three orders of magnitude for both 
simulations with Vine and Gadget, with both showing a 
remarkably similar evolution. The most notable difference to 
the simulations with direct magnetic field treatment shown 
in Fig. [5] and [6] is at the centre of the galaxies, where in 
the direct simulations the field amplification was strongest. 
With Euler potentials the magnetic field grows mostly in the 
spiral arms of the galaxy (see also Fig. [13]). 

Since the magnetic fields in our simulations are pas- 
sive, the density profiles (Figs. [2] and [3|) of the disc are the 
same for all runs. Thus, the different profiles of the magnetic 
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Figure 10. Numerical /i • | V • B|/|B| at t ~ 1.5 Gyr as a function 
of radius for simulations without applying the viscosity limiter 
using direct magnetic field treatment (blue) and using Euler po- 
tentials (black) in Gadget (solid lines) and Vine (dotted lines). 
Direct magnetic field simulations for which the viscosity limiter 
was applied are also shown (orange). Using direct magnetic field 
description, the numerical /i • V • B is highest at small radii, and 
much larger than in the Euler potentials formalism. 



field energy cannot be traced back to the density profiles. 
In fact, it is the numerical V • B which presumably causes 
the high amplification of the magnetic field at the centre in 
simulations with the direct magnetic field treatment. Fig. 
[TOl shows the radial profile of the numerical h • |V • B|/|B| 
at time t ^ 1.5 Gyr for simulations using direct magnetic 
field treatment (blue for simulations without applying the 
viscosity limiter and orange where the limiter was applied) 
and Euler potentials (black) performed using Gadget (solid 
lines) and Vine (dotted line). Utilising the direct magnetic 
field description, the numerical V • B is highest at small 
radii, and much larger than for the Euler potential formal- 
ism. As will be discussed in the following section, high V • B 
corresponds to high amplification of the magnetic field. 

Fig. II II shows the magnetic field vectors for the normal 
resolution Vine simulation utilising Euler potentials at the 
time t ^ 0.9 Gyr. This time the colours correspond to the gas 
density on a logarithmic scale from 0.3 • 10~^ to 2.3 • 10^ Mq 
pc""^, overplotted with the field vectors. The length / of the 
vectors is normalised to the initial value and displayed log- 
arithmically as / = 3 • \og(B/Bo), i.e. / = corresponds to 
B ^ Bo or smaller, I = 1 to B ^ 2 ■ Bo, I = 2 to B ^ 5 ■ Bq 
and I — 3 to B = 10 ■ Bo. The magnetic field lines follow 
the spiral structure of the gas. They have been amplified by 
contraction in regions of higher density and restructured by 
differential rotation of the galaxy. Their orientation is caused 
by the motion of the gas. These characteristics are very simi- 
lar to typical observations of magnetic fields in galactic discs 
(e.g. Fig.IU). 

Qualitatively, this behaviour is the same for all simu- 
lations using both codes. Only the central region in simu- 
lations using direct magnetic field treatment shows chaotic 
orientation of the magnetic field lines, indicating artificial 
amplification of the magnetic field due to high numerical 
V B. 
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Figure 12. Same as Fig. [G] this time the magnetic field is smoothed every 30 timesteps. The field morphology is similar to the 
morphology in Fig. [6] and [5] respectively, but the magnetic field values in the spiral arms are now more similar to the values in the Euler 
implementation. 
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Figure 11. Gas density (colour coded) and magnetic field vec- 
tors for the normal resolution simulation with Vine using the 
Euler potentials formalism at t ~ 1 Gyr, i.e. ~ 500 Myr after the 
inclusion of the magnetic field. The length of the vectors is nor- 
malised to the initial value and displayed logarithmically. I = 
corresponds to B ^ Bq or smaller and I = 3 to B = 10 • Bq. 



5 EVALUATION 

5.1 Magnetic field growth 

Figs. [5] dH and [7] JS]), respectively, reveal the differences in 
the magnetic field amplification for the direct magnetic field 
treatment and the Euler potentials formalism: Using the di- 



rect description, the amplification of the magnetic field en- 
ergy in the spiral arms is higher by at least two orders of 
magnitude, and at the centre even more than six orders of 
magnitude compared to the Euler potentials method. This 
difference is probably caused by the numerical V • B in these 
simulations fFig. I1Q|) , but possibly also by the fact that field 
winding is not traced beyond a certain evolutionary state 
in the Euler potentials formulation (see section [3 .1.6^ . Since 
the Euler potentials are free from physical divergence by 
construction (i.e. the divergence is zero to measurements 
errors), the numerical divergence in simulations using the 
Euler potentials is due to the SPH derivative approximation 
when calculating the magnetic field from the potentials (Eq. 
lisp . In this sense, the numerical divergence found in sim- 
ulations using Euler potentials reflects the ability of SPH 
operators to measure the gradient of a curl to zero. Thus, 
the fact that V • B is higher by approximately one order 
of magnitude in the disc (i.e. within 5 to 15 kpc) and 
by several orders of magnitude at the centre (Fig. [10]), pre- 
sumably causes the different magnetic field amplification in 
these simulations. This is the case at least in the disc re- 
gion, where the winding of the field is not strong enough to 
constrain the Euler potentials formulation. 

To get a better idea of the influence of numerical V-B on 
the amplification of the magnetic filed, we have performed 
simulations applying magnetic field smoothing, a tech- 
nique allowing for reduction of small scal e fluctuations and 
therefore also the numerical divergence ( Dolag Stasvszvnl 
I2OO8I ) . Within this method, t he magnetic f i eld is smoothed 
periodically as suggested by lB0rve et all (|200lh . Fig. [12] 
shows again the magnetic field energies and gas densities 
for a Gadget simulation starting from the same initial con- 
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Figure 13. Total magnetic field (Btot = yj Bl + Bl + B^) at 
t ^ 1.5 Gyr as a function of radius for the normal resolution 
Gadget (solid lines) and Vine (dotted lines) simulations with- 
out applying the viscosity limiter using Euler potentials (black) 
and the direct magnetic field description (blue). Direct magnetic 
field simulations for which the viscosity limiter was applied are 
also shown (orange). Two direct magnetic field implementations 
including field smoothing to reduce the numerical V • B contri- 
bution are shown in red and green. The simulations including 
smoothing have been run with a smoothing interval of 30 (red) 
and 5 timesteps (green), respectively. Increasing the frequency of 
smoothing tends to decrease the amplitude of the magnetic field 
between 3 and 10 kpc, but has relatively little effect for larger 
radii. 



ditions as before and without applying the viscosity lim- 
iter. This time, the magnetic field was smoothed every 30 
timesteps. Applying the smoothing scheme, the amplifica- 
tion of the magnetic field energy is reduced to approximately 
three orders of magnitude within the spiral arms, which is 
the same as the amplification seen in simulations using the 
Euler potentials, and it is also lowered towards the centre 
of the galaxy. The structure of the magnetic field is despite 
the smoothing still very similar to the other runs and again 
correlates well with the structure of the gas density, how- 
ever, the magnetic field energy is more concentrated within 
the spiral arms. 

Fig. 1131 shows the total magnetic field at t ~ 1.5 Gyr 
as a function of radius for the normal resolution Gadget 
(solid line) and Vine (dotted line) simulations using Euler 
potentials (black) and the direct magnetic field description 
without applying the viscosity limiter (blue) and with the 
limiter turned on (orange), respectively. The direct magnetic 
field simulations including field smoothing are shown in red 
and green. They have been performed with a smoothing in- 
terval of 30 (red) and 5 timesteps (green), respectively, and 
without applying the viscosity limiter. As discussed before, 
the most notable difference between simulations with direct 
magnetic field treatment and the Euler implementation is 
at the centre of the galaxies. There, the amplification in the 
direct simulations is much stronger than in the Euler simu- 
lations. This behaviour could be at least partly physical, as 
there are high radial velocities and strong in- and outflows of 
gas in the central region (fig. [17]), resulting according to Eqs. 
[71 and [8] in an amplification of the magnetic field. In addi- 
tion, also the azimuthal derivatives of the radial and toroidal 



velocity components are large at the very centre, which also 
could account for the violent amplification (see section 15. 3p . 
On the other hand, the second term of eq. [8] does not play 
an important role, since dv^/dr is large and therefore dQ/dr 
small in the central region (by reason of solid body rotation). 
However, the high V • B values at the centre make it diffi- 
cult to distinguish between physical growth and numerical 
errors. Since the Euler potentials are also unreliable in this 
region (see section [3. 1. 6 1) , it is not easy to decide which for- 
malism is the most capable in describing the physics in the 
centre of the galaxy correctly. This is also true for the sim- 
ulations including smoothing. Increasing the frequency of 
smoothing tends to decrease the amplitude of the magnetic 
field between 3 and 10 kpc, but has relatively little effect 
for larger radii. Interestingly, the large increase of B in the 
centre is never smoothed away, which could indicate, that 
this behaviour is actually partly physical. For the simula- 
tion which applies smoothing every 30 timesteps (red), the 
amplification of the field at r > 3 kpc is similar to the sim- 
ulations with Euler potentials. Applying smoothing every 5 
timesteps (green), the amplification is considerably weaker 
than in the Euler potentials simulations, indicating that by 
such strong smoothing essen tial physics is lost, in agreement 
with earlier findings by Dolag Stasvszv"nl (|2008[ ) . 

In the following, we only consider the disc region (from 
5 to 15 kpc), since the high numerical divergence in the 
centre makes it difficult to lower it to the value of the 
divergence seen in simulations with Euler potentials (i.e. 
h • |V • B|/|B| ^ 1), without smoothing the magnetic field 
structure too much. 

Fig. [TJ] shows the evolution of the total magnetic field 
{Btot = y^B|~+~B|~+~B|) within the disc with time for the 
different implementations. The colour coding is the same 
as in Fig. 1131 As before, for the simulation which applies 
smoothing every 30 timesteps (red), the amplification of the 
field is similar to the simulation with Euler potentials. How- 
ever, the performance of these simulation is not very con- 
vincing due to the "jumps" in the evolution caused by the 
artificial periodic smoothing. Applying smoothing every 5 
timesteps (green), the amplification is as discussed before 
lower than in the Euler potentials simulations. 

This behaviour can be understood by considering the 
corresponding numerical divergence of the magnetic field. 
Fig. [T5l shows /i-|V-B|/|B| as a function of time for all 
simulations. In all cases, the growth of h • |V • B|/|B| be- 
haves similar to the amplification of the total magnetic field, 
i.e. the higher the divergence, the stronger the amplification 
of the field. Though the numerical divergence in the simu- 
lation using Euler potentials (black) is higher than in the 
simulation with a smoothing interval of 5 timesteps (green) , 
its value does not directly correlate with the field growth. 
That is because the (defective) magnetic field itself is not 
used for calculating the magnetic field evolution within the 
Euler potential formalism as is the case for the direct mag- 
netic field description (compare eas. 1171 and I18|) . Using the 
smoothing scheme lowers the divergence (in case of smooth- 
ing every 5 timesteps even below the numerical divergence 
of the Euler potential formalism) and lowers also the field 
amplification, leading (if applied not too often) to an am- 
plification of the total field much more similar to that using 
the Euler potentials, which are free from physical divergence 
by construction. 
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Figure 14. Total magnetic field (Btot = ^ B'^ + + ^1) 
within the disc (between 5 and 15 kpc) as a function of time 
for different implementations. The colour coding is the same as 
in Fig. [13] Applying the smoothing scheme reduces the amplifi- 
cation of the magnetic field. 
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Figure 15. Total divergence of the magnetic field within the 
disc (between 5 and 15 kpc) as a function of time for different 
implementations. The colour coding is the same as in Fig. [T3l 
The higher the divergence, the stronger the amplification of the 
magnetic field. 



Interestingly , for simulation s applying the viscosity lim- 
iter suggested bv lBals ara (1995), the magnetic field amplifi- 
cation using the direct magnetic field description is in both 
codes much higher than without applying this limiter (or- 
ange lines in Figs. 1131 and [T4)) . The reason for this higher 
amplification is the higher velocity dispersion in these simu- 
lations. The viscosity limiter lowers the viscosity in regions 
of strong shear flows, thus suppressing velocity diffusion and 
leading to higher velocity gradients. Consistently, also the 
numerical divergence of the magnetic field is higher (and 
considerably higher than the "unavoidable" value of approx- 
imately one) in these simulations (orange lines in Figs. [TOl 
and [T5|) . Applying the viscosity limiter in simulations using 
Euler potentials, however, does not change the evolution 
of the magnetic field significantly (not shown). Therefore, 
again, it is probable that the higher numerical V • B terms 
lead via the induction equation (eq. [3]) to an enhanced mag- 
netic field growth. 

In summary, the field amplification in case of direct 
magnetic field description correlates with the non- vanishing, 
numerical V • B. The Euler potential formalism also has its 
shortcomings (like the non-uniqueness and the dependence 
of the magnetic field on two derivatives (Eq. I18|) leading 
to lower numerical accuracy). Thus there is a strong need 
for simulations with different V • B cleaning techniques and 
even higher resolution in order to be able to distinguish the 
best description for simulations of magnetic fields in galactic 
discs. 

However, since the physical divergence is zero in the 
case of the Euler potentials, we believe this method (for the 
time being) to be the best one for our studies of magnetic 
convection in disc galaxies. The following discussion there- 
fore concentrates on simulations using Euler potentials. 



5.2 Numerical resolution 

Fig. 1161 shows the total magnetic field as a function of time 
for different resolutions (see Table [2]) in simulations with 
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Figure 16. Resolution study: Total magnetic field (Bu 



+ + B|) as a function of time for the Gadget (solid 
line) and Vine (dotted line) simulations without applying the 
viscosity limiter using Euler potentials. The total numbers of par- 
ticles are 1.3 • 10^ (blue), 1.3 • 10^ (black) and 1.3 • lO'^ (red). 



Gadget (solid lines) and Vine (dashed lines) without ap- 
plying the viscosity limiter. 

One Gyr after its initialization the magnetic field has 
been amplified from 10~^ to approximately 9 • 10~^ G in 
the low resolution simulation (blue), whereas the final mag- 
netic field strength in the normal resolution simulation is 
slightly more than 1.5 times higher (1.5 • 10~^ G). The final 
magnetic field strength in the high resolution run is again 
approximately 1.5 times higher than in the normal resolu- 
tion run (i.e. ^ 2.5 • 10"^). The numerical /i • |V • B|/|B| 
values are of the same order for all resolutions (not shown). 
Thus, we have not yet reached numerical convergence in the 
magnetic field evolution. 
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5.3 Inspection of the induction equation 

By analyzing the velocity and magnetic field in our sim- 
ulation we can identify the single terms of the induction 
equation responsible for the behaviour of the magnetic 
field. Dropping all dependencies on the equations for the 
evolution of the radial and toroidal magnetic fields read 
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where we have labelled the single terms with numbers 
for easier reference. 

The radial and toroidal components of the velocity and 
the magnetic field and their corresponding derivatives are 
shown in Fig. [17] after approximately three half mass rota- 
tion periods after the onset of the magnetic field. The radial 
velocity (top left) is typically negative, leading to an effective 
gas inflow towards the centre of the galaxy. This negative 
radial velocity mirrors the angular momentum transport to 
large radii of the galaxy by spiral arm formation. The mean 
circular velocity (second row, left) is 210 km s~^ at large 
radii, and drops to zero towards the centre (see also Fig. 2]). 
The toroidal magnetic field (bottom left) is wound up by dif- 
ferential rotation, leading to a structure of altering positive 
and negative magnetic field values from centre to the edge of 
the galaxy. Consequently the derivatives with respect to cp 
(right panel) are smaller than the radial derivatives (middle 
panel), mirroring the approximate axial symmetry. However, 
since the terms of the induction equation depend always on 
a product between a derivative and a velocity or magnetic 
field component, one cannot a priori neglect the terms de- 
pending on azimuthal derivatives. 

In order to quantify the influence of the different terms 
1-10 during the simulation we calculated their values in 
cylindrical bins within the disc (5 to 15 kpc) and their mean 
value at different times. We have taken the negative values 
of each term in case of negative magnetic field to distin- 
guish between amplifying and attenuating terms. The result 
of this calculation is shown in Fig. [18] The upper plot shows 
the temporal evolution of the terms responsible for ampli- 
fication/attenuation of the radial magnetic field (terms 1 
to 5) and the lower of the toroidal magnetic field (terms 
6 to 10). Positive values imply amplification, and negative 
attenuation of the corresponding B- component. The non- 
axisymmetric terms are shown in red. 

Looking at Fig. llSi the most important term for the evo- 
lution of the radial magnetic field is term 5, i.e. 
Since the toroidal velocity dominates the velocity field, 
this term is most important although is comparatively 
small. This can be seen following the evolution of the circu- 
lar velocity and the radial magnetic field more closely: The 



radial magnetic field is strongest where the circular veloc- 
ity has its highest value, with a delay of roughly 40 Myr. 
All other terms lie in the same range and therefore com- 
pete with each other. Since their values are positive as well 
as negative, one should not expect a significant amplifica- 
tion on their account. This analysis shows, that even small 
deviations from axial symmetry are very important for the 
evolution of the magnetic field in spiral galaxies. 

One reaches the same conclusion looking at the terms 
of the evolution equation for B^. Except for the beginning 
of the simulation, the leading term here is clearly term 9, 
i.e. the only non-axisymmetric term in this equa- 
tion. Term 10, which was our candidate for the most im- 
portant term for axial symmetry, is only the second most 
important. Both terms depend on the toroidal velocity com- 
ponent, thus demonstrating the importance of the differen- 
tial rotation for the evolution of the toroidal component of 
the magnetic field. 

Neglecting all non-axisymmetric terms (plotted in red) 
one finds term 1 (—Br^) to be largely dominant over term 
term 4 {—Vr^§^), in agreement with the theory for the evo- 
lution of Br. Also the term responsible for the evolution of 
B^ is as expected: Term 10 { — ^Br) is the leading term 
and followed by term 6 {—B^ 



V dr 



). However, term 1 and 6 
are both of order 10"^"^ G Myr~^, thus not being able to 
account for any significant amplification of our initial mag- 
netic field, and term 10 can only amplify B^ effectively if 
Br is amplified. 

This behaviour is qualitatively the same also for runs 
with the direct implementation of the induction equation. 
We conclude that the non-axisymmetry of the system is the 
driving force for the observed field amplification in our sim- 
ulations. 



6 CONCLUSION AND OUTLOOK 

We have presented a set of self- consistent simulations of the 
evolution of magnetic fields in galactic discs performed with 
the AT-body codes Vine and Gadget. Hydrodynamics was 
treated using the SPH method. The evolution of magnetic 
fields within the framework of ideal MHD was followed by 
both a direct implementation of the induction equation and 
a formalism using Euler potentials. 

The presented set of simulations shows the importance 
of a sensible treatment of V • B when simulating magnetic 
fields in spiral galaxies. Since artificial magnetic monopoles 
can be responsible for unphysical amplification of the field, 
more studies of possibilities to avoid or inhibit numerical 
V • B terms are still needed. Although the description us- 
ing Euler potentials avoids (physical) magnetic monopoles 
by construction, the drawback in using them is that they 
lead to constraints on magnetic helicity. Since helicity fluxes 
can affect the dynamo process within a mean field the- 
ory ([Brandenburg Subramanianl l2005l l , Euler potentials 
would probably not be suitable for simulations including the 
a-effect. Furthermore, Euler potentials do not allow for all 
initial field configurations, since they are not necessarily sin- 
gle valued and in addition, their derivation can become quite 
complex. Nevertheless, using topologically simple initial con- 
ditions for the magnetic field, the Euler potential formalism 
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Figure 17. Radial and toroidal components of the velocity and the magnetic field and their derivatives at t ^ 900 Myr for the normal 



resolution Gadget simulation. From left to right and top to bottom: -Ur, ^ ^ • ^ '^^^ ~ ' 



dv^ dBr 1 dBr TD 



and 



1 _ dB^ 

r dtp 



seems to be the best tool to follow the ideal evolution of 
magnetic fields in simulations of spiral galaxies with SPH. 

A possible alternative to Euler potentials is the vec- 
tor potential A. The disadvantages are the need for a time 
integration of A when evolving magnetic fields and the oc- 
currence of second derivatives in the force equation when 
calculating magnetic forces (Fmag ocjxB(x(VxB)xB(x 
(V X (V X A)) X (V X A)), both leading to lower accuracy 



in the calculation. On the other hand, the advantages are a 
somewhat easier derivation of A for a given magnetic field 
and that there are no constraints on magnetic helicity using 
a vector potential. It would be definitely interesting to study 
the differences between simulations utilizing a vector poten- 
tial and the Euler description, although it could be hard 
to overcome the problems related to numerical intricacies 
within a SPH implementation of the vector potential. 
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plifying process should therefore account for four orders of 
magnitude of increase within few Gyrs in order to reach 
the observed values of ^ 10/iG. Since our simulations of a 
purely kinematic dynamo account at best for one order of 
magnitude, there is still need for a more complete scenario 
with additional subgrid physics. Such subgrid physics should 
comprise the a-effect due to turbulent gas motions below 
the resolution limit, estimated from local high-resolution 
MHD simulations and observations of turbulent motions 
in nearby galaxies. Hereby, potentially the most promising 
ansatz is the cosm ic ray driven dynamo (|Lesch Hanas j 
l200l [ Hanasz et al. Given the fact, that the presented 

simulations reveal the complete three dimensional velocity 
field to fully account for the large-scale structure of the mag- 
netic field, we believe that iV-body SPH together with sen- 
sible subgrid physics will be apt to test our understanding 
of the evolution of magnetic fields in spiral galaxies numer- 
ically. 




0.5 0.8 1.0 1.2 1.4 1.6 1.8 2.0 
t [Gyr] 

Figure 18. Values of the different terms of the induction equa- 
tion with time for the normal resolution Gadget simulation us- 
ing Euler potentials. Upper plot: temporal evolution of the terms 
responsible for amplification/attenuation of the radial magnetic 
field (terms 1 to 5). Lower plot: Evolution of terms 6-10, responsi- 
ble for the evolution of the toroidal magnetic field. Positive values 
imply amplification, and negative attenuation of the correspond- 
ing B-component. The non-axisymmetric terms are shown in red, 
the axisymmetric terms are shown in black. 



The analysis of the different terms of the induction 
equation applied to our simulations clearly show that the 
non-axisymmetry of the velocity and magnetic field cannot 
be ignored in any consideration of the kinematic dynamo. 
There are two main processes leading to angular momentum 
transport and hence non-axisymmetry in spiral galaxies: In- 
ternal driving due to spiral structure and bar formation (the 
former considered in the presented paper) and external driv- 
ing due to interaction with other galaxies. Simulations of 
interacting systems would therefore enrich our understand- 
ings even further on how large scale magnetic fields evolve 
due to large scale velocity fields. 

Our simulations show only a weak amplification of the 
initial magnetic field. Observations of spirals galaxies at high 
redshifts suggest that their magnetic field strengths were at 
least as strong as the magnetic fiel ds at the current epoc h 
within few Gyrs of the Big Bang ([Kronberg et all I2OO8I I. 
Assuming initial strengths of order Biqm — 10~^ G an am- 
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